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We numerically study the superfluid to Mott insulator transition for bosonic atoms in a one dimen- 
sional lattice by exploiting a recently developed simulation method for strongly correlated systems. 
We demonstrate this methods accuracy and applicability to Bose-Hubbard model calculations by 
comparison with exact results for small systems. By utilizing the efficient scaling of this algorithm 
we then concentrate on systems of comparable size to those studied in experiments and in the pres- 
ence of a magnetic trap. We investigate spatial correlations and fluctuations of the ground state 
as well as the nature and speed at which the superfluid component is built up when dynamically 
melting a Mott insulating state by ramping down the lattice potential. This is performed for slow 
ramping, where we find that the superfluid builds up on a time scale consistent with single-atom 
hopping and for rapid ramping where the buildup is much faster than can be explained by this 
simple mechanism. Our calculations are in remarkable agreement with the experimental results 
obtained by Greiner et al. [Nature (London) 415, 39 (2002)]. 

PACS numbers: 03.75.Lm 



I. INTRODUCTION 

Recent experiments on loading Bose-Einstein con- 
densates into an optical lattice have allowed for the 
creation and study of strongly correlated systems of 
atoms 0,HI1I3. In particular the superfluid (SF) to 
Mott insulating (MI) transition first observed in a semi- 
nal experiment by Greiner et al. [lj has received a lot of 
attention since it impressively demonstrated a clean real- 
ization of the Bose-Hubbard model (BHM) p| which has 
long been considered a toy model in condensed matter 
physics. Furthermore, in the ideal MI state each atom is 
localized to a lattice site corresponding to a commensu- 
rate filling of the optical lattice with zero-particle-number 
fluctuations. These properties make MI states attractive 
candidates for several applications,most notabl y q uan- 
tum memory, quantum computing 0, 0, H, 0, HH LUl LL2l j 
and quantum simulations of many-body quantum sys- 
tems 111 El- 

The BHM Hamiltonian describes atoms loaded into a 
sufficiently deep optical lattice 0, . It contains a ki- 
netic energy term, with matrix element J, describing the 
hopping of particles from one site to the next and an in- 
teraction term, with matrix element U, which accounts 
for the repulsion of two atoms occupying the same site. 
The ratio U/J increases with the depth of the optical 
lattice and can be varied over several orders of magni- 
tude by tuning the optical lattice parameters ^5|. In 
particular by changing the intensity of the laser beams 
creating the optical lattice it is possible to vary J and U 
on time scales much smaller than the decoherence time 
of the system. This opens up the possibility of directly 
studying the dynamics of the BHM during the quantum 
phase transition at temperature T = According 
to mean-field (MF) theory this phase transition occurs 
at u c = U/zJ ps 5.8, where z is the number of nearest- 
neighbor sites in the lattice [1 0, and is easily ac- 
cessible in an optical lattice. 



In [l| the dynamics of atoms in a three dimensional 
(3D) optical lattice (z = 6) was studied while more 
recently optical lattice setups where the motion of the 
atoms was restricted to ID {z = 2) @ were investigated. 
These experiments revealed some striking properties of 
the quantum phase transition. In particular a feature 
which is yet to be fully understood is the time scale over 
which coherence is built up throughout the atomic system 
when going from the MI to the SF limit [lj . Indeed it can- 
not be easily explained using MF theory and numerical 
studies of this dynamical effect were, until now, limited 
to small systems of approximately ten atoms. Recently, 
however, it has been shown that quantum computations 
on ID systems of qubits which do not give rise to strong 
entanglement can be efficiently simulated on a classical 
computer via the so called time-evolving block decima- 
tion (TEBD) algorithm 20]. An immediate application 
of this discovery is to the simulation of the time evolution 
of many-body ID quantum systems which are governed 
by a nearest-neighbor Hamiltonian [2l| . The BHM is one 
of many important model Hamiltonians which fall into 
this class [23. The simulation method is efficient for all 
such ID model Hamiltonians due to a universal property 
of ID systems that their ground state and lowest-lying 
excitations tend to contain only a small amount of en- 
tanglement |2l| . 

In this paper we restrict our attention to the ID BHM 
with our physical motivation being to study the the na- 
ture and speed at which the superfluid component is 
built up as the system is dynamically driven through the 
SF-MI transition. By exploiting the efficient scaling of 
the TEBD algorithm with the size of the system we are 
able to investigate this phenomenon for setups which are 
of comparable size to those studied in experiments 0. 
First, in Sec. [H] we introduce the ID BHM for describ- 
ing atoms in optical lattices and briefly introduce the 
TEBD algorithm as used in this paper. In Sec. IIIII we 
then demonstrate the applicability of the TEBD to the 
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BHM by comparison with exact numerical calculations 
for small systems. This is then followed by an investiga- 
tion of SF and MI ground states of larger lattice setups 
concentrating on their spatial correlations and occupa- 
tion number fluctuations together with a comparison to 
MF results. We then study the dynamics of the MI to SF 
transition in Sec. II VI when changing the lattice depth on 
two different time scales. Most notably for rapid MI to 
SF ramping we find that the width of the central inter- 
ference peak, as observed after releasing the atoms from 
the lattice, shrinks with an increasing total ramping time 
with the same functional dependence found in 0. This 
result is discussed in Sec. IIVI B. Finally, we summarize 
our results in Sec. E] 



II. MODEL AND NUMERICAL METHOD 

In this section we introduce the BHM describing 
bosonic atoms in an optical lattice where the motion is 
restricted to ID and give a short overview of the numer- 
ical method used in our simulations. 



A. Model 

By confining an ultracold bosonic gas in a 3D optical 
lattice with a large depth in the two orthogonal directions 
y and z it is possible to create an array of effective ID 
systems in the x direction 0, IToL |23| . The dynamics 
of these systems is governed by the external trapping 
and the optical lattice potential along the x axis. The 
optical lattice then has a depth Vb proportional to the 
laser intensity and a lattice period a = A/2, where A 
is the wavelength of the laser light. The Hamiltonian 
describing each ID system reduces to the ID BHM (for 
details see Appendix [XJ (taking H = 1 throughout): 

H = - J(6j 7t 6 m+ i + H.c.) - [i m b\j> m + —blnblnbmbm, 

m ~ 

(1) 

where the operators b m (b^) are bosonic destruction (cre- 
ation) operators for a bosonic particle in site to, centered 
at x m — ma, obeying the standard canonical commu- 
tation relations. The grand canonical Hamiltonian then 
has n m = ii — Vr(x m ) as the local chemical potential 
for site to, where Vr is the external trapping potential. 
The parameters U and J can be determined in terms 
of the Wannier functions w(x) as shown in Appendix IaI 
and under the assumptions outlined are independent of 
the lattice site to |15|. Their ratio can be varied over a 
wide range by dynamically changing the depth Vo of the 
optical lattice. For all the systems considered here we 
take the wavelength of the light used to form the optical 
lattice as A = 826 nm and the atomic species trapped as 
87 Rb, where a s = 5.1 nm. 



B. Numerical method 

In this paper we exploit the recently devised TEBD 
simulation algorithm [53, which allows the dynam- 
ics of ID systems with nearest-neighbor interactions, 
such as the BHM, to be computed accurately and ef- 
ficiently. The TEBD algorithm has been shown to be 
closely related to the density matrix renormalization 
group (DMRG) [H |2^. Over the past decade the 
DMRG has provided enormous insight into the static and 
dynamic equilibrium properties of ID systems. Although 
originally devised as a ground-state method, it has been 
extended to yield accurate low-energy spectra |26j |. and 
also to calculations of the real time evolution of ID sys- 
tems [23 which is of particular importance here. The ap- 
proach of H3 is to take the DMRG ground state |V>(°)) 
obtained for the initial Hamiltonian and use it to define a 
decimation of the Hilbert space in which the Schrodinger 
equation is numerically integrated. The key assumption, 
and most severe approximation, within this scheme is 
that this static subspace defined by |?/>(0)) is adequate 
to approximate \ip(t)) with reasonable accuracy for all 
times. In general this will only be true for short periods 
of time. Novel methods have been devised j2^ which can 
maintain the accuracy over longer periods by 'targeting' 
other states in addition to the ground state, but in do- 
ing so the efficiency of the computation is significantly 
reduced [2j| . In contrast the TEBD algorithm can main- 
tain typical DMRG accuracies while remaining efficient. 
Despite their differing origins it has recently been shown 
that TEBD and DMRG algorithms share some crucial 
conceptual and formal similarities [22I I29I ] . Indeed both 
methods search for an approximation to the true wave 
function within a restricted class of wave functions which 
are described by matrix product states [30l l3l| , and do 
so with identical decomposition and truncation proce- 
dures. The essential difference, which we shall emphasize 
shortly, is that the TEBD algorithm updates the matrix 
product decomposition directly and in such a way that 
the resulting decimated subspace in which the time evo- 
lution is computed is optimally adapted at each step |22| . 

Here we briefly outline the essential features of the 
TEBD algorithm, with specific attention to its applica- 
tion to the BHM. Let us consider a ID BHM composed 
of M sites. An arbitrary state of this system can be 
expanded in the Fock basis 

00 00 

IV-') = ' ' ' c »i-»m I rai> ■ ■ ■ ,n M ) , (2) 

711— riM—0 

where \n m ) denotes the Fock state of n m particles in site 
to. For the purpose of simulating this system the number 
of Fock basis states per lattice site must be cut off to 
some upper limit n max - In all the numerical calculations 
we performed n max = 5. This is sufficient to avoid any 
cut off effects in the bosonic occupation, as long as only 
small filling factors of the lattice are used and the on-site 
interaction energy U is sufficiently large compared to the 
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FIG. 1: (a) The sequence of contiguous partitions of the sys- 
tem in which the SD are computed. The coefficients and 
states from these SD are then used to form the F and A ten- 
sors, (b) A depiction of the F tensors associated to lattice 
sites and A tensors associated to bonds between those sites. 



hopping energy J. 

Now suppose we split the system into two contiguous 
parts A m composed of the first m sites and B m composed 
of the last M — m sites. We can think of this partitioning 
as cutting the m th bond situated between sites m and 
m+ 1. For any state \tp) a Schmidt decomposition (SD) 
can be performed which renders the state in the form 



i^> = EAL m] K m >kf m > 

a = l 



(3) 



where \m is the Schmidt rank of the SD, Aq are the 
Schmidt coefficients, and \4> c a ), with c £ {A m , B m }, are 
the corresponding Schmidt states of the respective sub- 
systems. The Schmidt rank Xm is a useful measure of 
the entanglement between the two subsystems A m and 
B m H3- Given any state \ip) a set of (M — 1) SD can 
be performed according to a sequence of such partitions 
of the system with m £ {1- ■ ■ M — 1}, as depicted in 
Fig. da). 

Using the and states \4> s a ) for each subsystem ob- 
tained from these SD it is possible |23| to construct a set 
of F and A tensors which are equivalent to a matrix prod- 
uct decomposition of the expansion coefficients c ni ... nM 



of in the fixed Fock basis 

c n 1 ---n M = y ] 



Specifically one finds 



WiWAWri 2 ]" 2 aI 2 1 • • • \iM-i] v mn M 



( 4 ) 

where n m is the occupation number of site m, and a m 
are the Schmidt indices of the m th partition, each of 
which sums from 1 to its respective Schmidt rank Xm- 
With reference to Fig. ^b) we note that each A^ is 
labeled by the bond between sites m and m + 1, along 
with the corresponding Schmidt index a m , whereas each 
r|^"7a m is labeled by a site m which resides between the 
two bonds m — 1 and m, and so also possess the Schmidt 
indices a m _i and a m of these bonds. 

Under the circumstances described the expansion, 
Eq. , is exact and as such the number of parame- 
ters stored could grow exponentially with the size of the 
system. However, it is a general feature of ID systems 
with nearest-neighbor interactions that the entanglement 
within their ground state and low-lying excitations de- 
pends weakly on the size of the system |2l|. Indeed it 
can be shown that the entanglement of a block of size £ 
with the rest of the system remains finite as t — > oo in ID 
systems or at worst grows logarithmically with i at crit- 
icality £ 22j . Consequently the entanglement between the 
blocks of any of the (M — 1) SD illustrated in Fig.QJa) 
can be saturated by some fixed Schmidt rank, which for 
the systems we consider is typically small. It is this fact 
that accounts for the success of DMRG in ID systems. 
Similarly within TEBD it allows the maximum possible 
Schmidt rank used in the matrix product decomposition, 
Eq. (0J, to be fixed to some value %, thereby truncating it 
to the most significant contributions. For an appropriate 
choice of \ this approximation will be accurate with the 
error proportional to the sum of the discarded eigenval- 
ues in the SD |2jJ . This clear interpretation of the cen- 
tral numerical parameter x within TEBD is very useful. 
Once a value of x is found to saturate the entanglement 
of the ground state and low-lying excitations of a system 
then this a direct measure of the role of entanglement in 
the dynamics of the system. In total the scaling in the 
number of parameters within the expansion, Eq. is 
quadratic in x an d linear in the size of the system M 
and in rt ma x- So upon fixing x an d thus preventing its 
possible exponential dependence on M, the description 
becomes efficient. As with DMRG this decomposition of 
a state generates, for all practical purposes, an optimal 
X x x matrix product state |22|. A noteworthy limit of 
this is the approximation is where x = 1? which forces 
the description of the system to be of product form with 
respects to all sites. Using the TEBD algorithm under 
this severe restriction is in fact equivalent to MF theory 
and the Gutzwiller ansatz [1 [l7L 111 l3l l3^| . 

Another crucial advantage of the TEBD algorithm 
is that once a state is expressed in the matrix prod- 
uct form, Eq. one- and two-site unitary transfor- 
mations can be applied directly and exactly to the sys- 
tem such that the resulting state can be efficiently re- 
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turned to a matrix product form |22(. Indeed, given 
a partitioning of the system into a two-block-two-site 
configuration [1 • • • m — 1] [m m + 1] [m + 2 • • ■ M] , the 
application of a two-site unitary to sites m and m + 1 
only requires updates to be performed on the tensors lo- 
cal to those sites-namely, r^"™ Qm A^r[7^ 1 l™T +1 • The 
major computational effort of this update is limited to 
the rediagonalization of the reduced density matrix of 
one of the adjacent site and block subsystems, such 
as sites [m + l][m + 2 • • ■ Ml , which is of dimension 
(x n max) x (x^max) at most |20j . The crucial feature here 
is that a DMRG-style truncation to only the most rele- 
vant eigenstates of this reduced density matrix occurs in 
an optimal way at each application of a two-site unitary. 
This is in contrast with time-dependent DMRG methods 
where the basis states which make up the matrix product 
decomposition are fixed at the start p2l ? ] . The number 
of basic operations required to perform this update scales 
as 0( X 3 ) M- 

To compute the action of the time-evolution unitary 
exp (—iHdt), for a time step 5t, we first make the obser- 
vation that for Hamiltonians with nearest-neighbor in- 
teractions, which are composed of two-site operators at 
most, terms can be separated into a sum of those involv- 
ing odd sites, F, and those involving even sites, G: 

n odd 

G = G njn _|_i, (6) 

n even 

H = F + G. (7) 

Given that no terms within F involve the same lattice 
sites they all commute amongst themselves. Thus the 
action of exp (—iF8t) can be computed exactly as 

e -iFSt = Yl e -iF„, n+1 « (8) 
n odd 

Since each term in this product is a two-site unitary, they 
can be applied individually to the state with the method 
detailed in [23], and the same is also true for G. The 
complications in computing the time evolution arise from 
the fact that F and G do not in general commute, and 
hence we approximate the unitary time-evolution opera- 
tor exp(i(_F + G)St) using a Trotter expansion. Ignoring 
their noncommutativity would constitute a first-order ex- 
pansion. If we define 

S 2 (F, G, y) = e -iFy/2 e -iGy e -iFy/2^ (g) 

then the second-order expansion follows when y = St. 
For the numerical simulations performed in this paper 
the fourth-order expansion |34| was used, which has the 
form 

5 

e -i{F+G)St = -Q s ^ Gj qi§ ^ + Q ( St 5^ (1 0) 
1=1 



where the parameters qi are defined as 

qi = q2 = q± = 95 = q = 7^—^77^ 93 = 1- 4g. (11) 

A detailed analysis of the errors and computational cost 
of TEBD is given in |2j|, where it is shown that the 
Trotter error propagates quadratically with the simulated 
time and so the accuracy of the method can be main- 
tained for long periods with appropriate choices of the 
parameters. 

The pure TEBD implementation we employ here can 
be improved further by combining the advantageous fea- 
tures of TEBD outlined with the well-established op- 
timizations of DMRG such as good quantum numbers 
and White's 'state prediction' method. In doing so 
an adaptive time-dependent DMRG algorithm is ob- 
tained |22l . |29| illustrating the extremely close relation- 
ship between these two methods. Finally we note the 
very recent advances in generalizing TEBD and DMRG 
to describe mixed state dynamics and generic master 
equation evolution of ID systems with nearest-neighbor 
coupling [3^, |3|| . This opens up the possibility of sim- 
ulating finite temperature effects, decoherence and dissi- 
pation. 

III. GROUND STATES OF THE BHM 

We first investigate the ground state of the BHM and 
compare the numerical results with the exact ground 
states for a small homogeneous system. Then we consider 
a larger system in the presence of a shallow magnetic trap 
Vt superimposed on the lattice and compare the results 
to those predicted by MF theory. In all cases the numeri- 
cal ground state was computed with the TEBD algorithm 
using continuous imaginary time evolution from a simple 
product state, as detailed in |21| . 

A. Comparison of exact and simulated ground 
states 

To investigate the accuracy of the numerical simulation 
and its applicability to the BHM we first consider a small 
system in which an exact solution can be found readily. 
Specifically we use an optical lattice composed of M = 7 
sites, a trapping potential of Vt = with box boundary 
conditions, and a total number of particles N — 7. The 
ground state is then calculated numerically and exactly 
for U/2J = 2, 6 and 20, corresponding to the SF, in- 
termediate, and MI regimes, respectively. The numerical 
simulation was performed for x — 3, 5 and 7 in each case. 

The one-particle density matrices p m ,n — (bln^n) ob- 
tained for each regime for the numerical and exact cal- 
culations are visually indistinguishable in all cases. In 
order to highlight the extent of the agreement we present 
a number of other plots. Specifically in the SF regime 
the comparisons of the spatial correlation of the central 
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site |p4.4+d| as a function of the distance d are shown in 
Fig-I2fa) and Fig.|2Jb) between the exact and numerical 
calculations using x = 3 and x = 5. Identical compar- 
isons of the standard deviation of the site occupation 
<r(fim,m) = «JVm) ~ {N m ) 2 ) 1/2 , where N m = b\ n b m , and 
the spectrum of the one-particle density matrices e 7 [nor- 
malized with tr(p) = N] are shown in Fig.Gfc), Fig. Etd) 
and Fig. |2|e),Fig. |2Jf) respectively. These results show 
that although there is qualitative agreement between ex- 
act and x = 3 calculation, almost all expectation val- 
ues have a maximum deviation from the exact calcula- 
tion improved by an order of magnitude with % = 5. 
As expected for a SF ground state the one-particle den- 
sity matrix spectrum in Fig. Eff) is dominated by one 
eigenvalue of order N. However, given that the lattice 
still has a nonzero depth this state deviates from that 
of a pure SF, where IV'Sf) oc (J2 m ^ti) N l vac )j since the 
sum of the remaining eigenvalues (in descending order) is 
X^7=2 e 7 ~ 2.5 and so represents a significant quantum 
depletion of the SF. 

For the intermediate and MI regime a similar factor 
of improvement can be obtained; however in this case 
the x — 3 calculation already yields excellent agreement 
with the exact calculation. Specifically we find the infi- 
delity between the numerical and exact many-body state 
is 1 - F = 1 - | <^o) I ~ !0" 5 , where |0 O > (\%)) is 
the numerical (exact) ground state and the temperature 
corresponding to the difference in their ground-state en- 
ergy as AT w 10 _2 nK. The comparisons of the spa- 
tial correlation, site occupancy standard deviation, and 
one-particle density matrix spectrum for U/2J = 6 and 
U/2J = 20, with x — 3, are shown in Fig. [31 

These plots, along with those where x = 5 for the 
SF regime, illustrate the onset of increasing MI charac- 
teristics in the ground states. In particular the rapidly 
decreasing spatial correlations and fluctuations in occu- 
pancy, as well as the change in the spectrum from being 
dominated by one single-particle state to having seven 
almost equally occupied orbitals. These indicate that 
the MI ground state obtained is a very close approxima- 
tion to that of a pure MI, where \ipMi) oc J7 m |vac), 
representing commensurate filling of the lattice. However 
given that the lattice is not infinitely deep deviations with 
this pure MI state exist and are evident from the persis- 
tence of small off-diagonal correlations visible at d = 1 
in Fig. Ofb) and the spread of the spectrum about unity 
in Fig. Elf). 

As expected we find that the agreement between the 
numerical and exact calculations for a given value of x 
improves with increasing U/2J in line with the decrease 
in off-diagonal correlations. In all cases the x = 7 results 
gave excellent agreement with the exact calculation. The 
worst case being in the SF regime where an infidelity of 
1 — F w 10 -4 , and a deviation in ground-state energy of 
AT ps 10~ 2 nK was obtained. 

We note here that the exact calculation was performed 
with the canonical Hamiltonian with N = 7, whereas the 
numerical simulations used the grand canonical Hamilto- 



nian. With an appropriate choice of the chemical poten- 
tial \x the average particle number can be fixed to N = 7, 
enabling the comparisons above. Indeed for the calcula- 
tions performed the worst absolute value of the projection 
of the simulated state outside the TV = 7 Fock subspace 
was (00 1 Vnj=7 IV'o) ~ 10~ 13 . Hence our results confirm 
the agreement of these approaches for small systems, and 
we assume the agreement holds for the larger systems. 
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FIG. 2: Comparisons of the numerical (o) and exact (*) cal- 
culations with U/2.J — 2 for spatial correlations [p4,4+d| with 
the central site m = 4 obtained for (a) % — 3 and (b) \ — 5, 
the standard deviation of the site occupation a(p m ,m) ob- 
tained for (c) x — 3 and (d) \ — 5, and the spectrum e 7 of 
the one-particle density matrix obtained for (e) \ — 3 and (f) 
X = 5. The dashed and dotted curves shown are to guide the 
eye. 



B. MI and SF states with a superimposed 
magnetic trap 

To consider systems closer to those studied in experi- 
ments we use a lattice with M = 49 sites and made 
the system inhomogeneous by superimposing a harmonic 
trap potential Vr(x m ) — mAUJ 2 x' 2 n /2, where u> is the 
trapping frequency and mj is the mass of an atom. As 
with the smaller system the ground states for the SF, 
intermediate and MI regimes were calculated. The lat- 
tice was loaded with a total number of particles N = 40 
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FIG. 3: Comparisons of the numerical (o) and exact (*) cal- 
culations, where the numerics were all performed with \ = 3, 
for spatial correlations |/04,4+d| with the central site m — 4 
obtained with (a) U/2J = 6 and (b) U/2J = 20, the stan- 
dard deviation of the site occupation a(p m ,m) obtained with 
(c) U/2J = 6 and (d) U/2J — 20, and the spectrum e 7 of the 
one-particle density matrix obtained with (e) U/2J = 6 and 
(f ) U /2J = 20. Note the differing scales and that the dashed 
and dotted curves are shown only to guide the eye. 
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by choosing an appropriate chemical potential /i in each 
regime. For all cases a trapping frequency of oj = 2ir x 70 
Hz was used and found to be sufficient in eliminating any 
occupation at the boundaries of the system. 

The inhomogeneity caused by a spatially varying con- 
fining potential can result in the coexistence of spatially 
separated SF and MI regions. Such properties have 
been confirmed experimentally Q, 0] and have received 
intense theoretical study with numerical calculations. 
In particular through Gutzwillcr ansatz and MF the- 
ory [H H3 and Quantum Monte Carlo (QMC) [M H3 
and DMRG [4(j simulations in ID, as well as calculations 
using the QMC 0, ^2 an d numerical renormalization 
group 01 methods for 2D and 3D systems. Here we 
explore the SF-MI coexistence features of BHM ground 
states in order to confirm the physical picture arising 
from our numerical calculation for the large inhomoge- 
neous system. Specifically we make comparisons of the 
mean site occupancy and its standard deviation between 
the numerical results and those obtained by MF theory 
for each regime (see Appendix [B] for details of the MF 
calculation) [33ll3^| . 

The one-particle density matrix of the resulting SF 
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FIG. 4: The absolute value |p m ,n| of the one-particle density 
matrix as a function of site indices m, n for (a) U/2J = 2, (b) 
U/2J = 6, and (c) U/2J = 20. 



ground state is shown in Fig. Eta). Important features 
of this state are outlined in Fig. [5|a)-Fig. p%d). where it 
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FIG. 5: Specific plots for the three regimes. For U/2J = 2 
there is (a) the site occupancy |p m , m |, (b) the standard devi- 
ation of the site occupancy o"(p m ,m), both with the MF calcu- 
lation shown as the dashed curve, (c) the spatial correlations 
from the central site m = 25, and (d) the spectrum e 7 of 
the one-particle density matrix, showing only the 40 nonzero 
eigenvalues. The same plots are presented for U/2 J = 6 in 
(e)-(h) and for U/2J = 20 in (i)-(l). 



can be seen that the system is entirely SF. The MF re- 



sults for the site occupancy and its standard deviation, 
are also shown in Fig. 0a) and Fig. 0b), and as expected 
there is good agreement between the curves in both cases. 
For the intermediate regime, whose one-particle density 
matrix is shown in Fig. |4Tb). we see a system possess- 
ing alterna ting regions of coherent SF and incoherent MI 
phases [3^, |40| . The pattern of these regions starts with 
a central SF region with a mean occupancy exceeding 
unity, which then becomes a singly occupied MI, a SF 
region with mean occupancy less than unity, and finally 
the vacuum MI. In Fig.0e) and Fig.0f) we see that the 
MI region in Fig. 0e) coincides with suppressed fluctu- 
ations in occupancy shown in Fig.0f). The MF curves 
also plotted show general agreement with these phase 
identifications; however the MF curve in Fig. 0f) pre- 
dicts a significantly greater suppression of the particle 
number fluctuations for the MI regions than the numer- 
ical results. Such deviations are consistent with the fact 
that MF theory predicts a sharp and well-pronounced 
SF-MI phase transition [3^, • While these prediction 
are known to be accurate for infinite homogeneous sys- 
tems, for small inhomogeneous systems we see that the 
role of correlations is important and that the transition 
between the SF and MI regions is not established with 
such definitencss. 

Given that the ground state for the intermediate 
regime exhibits significantly sized SF regions which are 
separated by a MI, we investigated whether any corre- 
lations were present between these regions. Evidence 
of such correlations would be the presence of elongated 
peaks in the one-particle density matrix at the off- 
diagonal locations corresponding to the intersection of 
the rows m and columns n of the SF regions. To within 
the accuracy obtained with x = 5 no such correlation 
peaks were found in the one-particle density matrix of 
the ground state. This was confirmed with continuous 
imaginary time evolution not only for a product initial 
state used conventionally, but also with an initial state 
which already contained significant correlations. Our re- 
sults suggest that if such SF correlations do exist within 
the ground state, then they are extremely fragile. Despite 
this we shall see shortly that such correlations do occur 
readily in dynamical situations which cross the MI-SF 
transition, where the system does not necessarily remain 
in the ground state. 

Lastly, the one-particle density matrix for the MI 
regime is shown in Fig. Hfc). It is clear from this and 
the corresponding plots in Fig. 0i) and Fig. 0j) that 
the system is almost completely in the singly occupied 
MI phase, aside from the small SF regions at the far 
extremes before the vacuum. Their presence is typified 
by the two peaks in the occupancy standard deviation 
shown in Fig.0j). The MF calculation in this case gives 
the same identification of the regions, except for the very 
center of the trap, where a small SF region is predicted 
to exist, as evidenced by the central peak in the MF oc- 
cupancy standard deviation curve of Fig.0j). Again MF 
theory predicts a much greater suppression of occupancy 
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fluctuations in the MI regions of the system than the 
numerical calculation. As with the smaller system the 
plots of Fig. for the three ground states illustrate the 
transition from predominately SF to MI characteristics. 

In all cases the ground-state calculations were per- 
formed with x = 5. To ensure convergence the cal- 
culations were repeated for \ — 1 . The largest devia- 
tion between these two calculations was found in the SF 
regime where the estimated ground-state energy differed 
by a temperature AT X= 5^7 w 0.2 nK. We made a similar 
comparison between the x = 5 results and those of MF 
theory, which are equivalent to \ = 1, where the largest 
deviation was found to be AT x= i^$ w 3 nK. Given the 
larger occupancy of the system the calculations were also 
repeated with larger values of n max , confirming that no 
cut off effects were encountered. 



IV. DYNAMICS OF THE BHM 

The most novel feature of the TEBD algorithm is its 
capacity to efficiently simulate the dynamics of ID sys- 
tems which are inaccessible to exact calculation. Here we 
consider dynamics which are generated by varying the 
optical lattice depth Vb(i) according to some ramping 
profile in time. The results of this is that the parameters 
J(t) and U(t) in the BHM Hamiltonian of Eq. becom- 
ing time dependent. By appropriately choosing the range 
of values covered by the optical depth Vo(t) it is possi- 
ble to dynamically drive the system through the SF-MI 
transition. We shall consider such dynamics occurring 
on two different time scales: first via a slow and smooth 
profile, and second as fast linear ramping. Our objective 
being in both cases to observe the nature and speed in 
which coherence is reestablished within the system. 



A. Slow dynamics 

1. Profile of the slow dynamics 

For the slow dynamical profile a smoothed 'box' func- 
tion was used for the depth Vo(t). Such a profile for the 
dynamics has been considered before for small systems 
H3- It has the form 

V°V = Vsr+X 1 + ^-J% J/ti , (12) 

where t c , t w , and t s are time parameters specifying the 
center, width, and step size of the profile respectively, 
while M = 1 + e a is the scaling factor required to 
ensure that the depth varies from Vsf to Vmi- The lat- 
tice depths Vsf and Vmi were chosen to be the depths 
equivalent to U/2J — 2 and U/2J = 20, respectively, in 
correspondence with the parameters used in the previ- 
ous section for the SF and MI regime. The exact shape 
of the profile for U/2J resulting from the ramping of 



Vo(t) chosen is shown in Fig. Efa). Large time parame- 
ters have been used in order to keep the time evolution of 
the system sufficiently adiabatic and to prevent excessive 
excitations. 



2. Slow dynamics of the small system: M — 7 

First we consider the slow dynamics applied to the 
small system. This provides the opportunity to solve the 
BHM dynamics both numerically and exactly, allowing a 
direct comparison of the accuracy and a demonstration 
of the applicability of the algorithm to the dynamics of 
the BHM. The system was initially prepared in the SF 
ground state computed earlier in Sec. IIIII A, and using 
X = 5 for the numerical calculations. The time evolution 
was then performed for a total time ttot = 2t c , with time 
t running over the interval [0,t to t]- The spectrum e 7 of 
the one-particle density matrix p m>rl (t) is plotted as a 
function of time in Fig. HJb). For times t where U/2J < 
u c the spectrum is, as expected, dominated by one large 
eigenvalue whose value is of order of the number of atoms. 
As U /2J crosses the MF critical value u c the eigenvalues 
are found to converge around the region of 1. Indeed the 
state of the system given by the numerical calculation 
at the time t = t c in the dynamics, which corresponds 
to U/2J = 20, is found to have an infidelity with the 
numerical MI ground state computed earlier in Sec. IIIII A 
as 1 — F < 10~ 4 . This confirms that for a small and 
homogeneous system the ramping is sufficiently adiabatic 
to ensure that the system has entered the MI regime as 
the ground state and the one-particle density matrix is 
diagonal. 

With decreasing optical depth and, in turn, decreasing 
U/2J, the SF ground state is restored when U/2J = 2 
is reached again at t = 2t c . The infidelity between the 
initial numerical SF ground state and the final numerical 
SF state was found to be 1 - F < 10~ 3 . In Fig. EJc) 
the behavior of the fluctuations in the site occupancy 
is as expected; namely, the standard deviation in site 
occupancy is suppressed with increasing lattice depth, 
and restored with its subsequent decrease. 

To test the accuracy of the TEBD algorithm a number 
of comparisons to the exact calculation were made. The 
simplest of these was the maximum relative deviation be- 
tween the exact and numerical results for the one-particle 
density matrix spectrum Ae = max 7 (|l — e T /e' |), where 
e 7 and e 7 are the numerical and exact results, respec- 
tively. The time profile Ae is plotted in Fig. Eld). It 
is found that over the whole time evolution the relative 
deviation is at most Ae ~ 10 _1 . A similar relative devi- 
ation can be defined for the standard deviation of the oc- 
cupancy as Act = max m ( 1 1 — o~ m / a' m | ) , where it is found 
that Act w 10 -2 at most during the time evolution. 

The most conclusive comparison, however, is the in- 
fidelity 1 — F(t) of the exact and numerical many-body 
states over the time evolution, shown in Fig. GJe). It is 
clear from this that the infidelity is bounded as 1 — F(t) < 
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4 x 10 -3 over the whole evolution. The shape of the infi- 
delity profile also gives important information about the 
TEBD method. Namely, it fits the general observation 
made in Sec. lIIIl A for the ground-state calculations that 
for fixed numerical parameters x an d "max the simula- 
tion is more accurate in the intermediate and MI regimes 
than the SF regime. This behavior is precisely exhibited 
in the time dependence of the above comparisons where 
significant reductions in the deviations are seen when the 
system enters the MI regime. 




FIG. 6: Slow dynamics of the small system: (a) the result- 
ing ramping profile of the parameter U/2J with time, where 
the time parameters for the Vo(i) profile are t c = 60 ms, 
t w — 24 ms, and t s = 18 ms, (b) the spectrum of the one- 
particle density matrix e 7 with time obtained from the nu- 
merical calculation, (c) the standard deviation of the site oc- 
cupancy a(p m ,m) with time obtained from the numerical cal- 
culation, (d) the maximum deviation between the numerical 
and exact spectrum Ae with time, and (e) the infidelity 1 — F 
of the numerical many-body state compared to the exact state 
with time. 



3. Slow dynamics of the larger system: M = 49 

The larger system possesses many of the essential char- 
acteristics of current experimental implementations of 
the SF-MI transition 0,0 In particular these include 
the inhomogeneous nature of the system, caused by a 
trapping potential, and the larger number of both lattice 
sites, and atoms, as compared with the smaller system. 
With a linear size of M = 49 sites the system consid- 
ered is on the same scale as experiments already per- 
formed |lj,|2|- ^ ne ma j° r difference is that the mean oc- 
cupancy at the center of our system is (n c ) « 1.5, roughly 
half that of most experiments, where it is usual to have 
(n c ) as 2.5. While the mean occupancy undoubtedly has 
an important influence in the dynamics, the system sim- 
ulated here is sufficiently close that it can demonstrate 
much of the important physics. 

The slow ramping profile, Fig. IHIa), was performed 
identically on the larger inhomogeneous system using 
the ground state Fig. Ufa) computed earlier as the ini- 
tial state and for a total time ftot = 3i c . The resulting 
spectrum e 7 of the one-particle density matrix is plotted 
as a function of time in Fig. General features of this 
spectrum follow from the smaller homogeneous system; 
namely, the trend for the eigenvalues to decrease as the 
bottom of the ramping is reached. While most eigenval- 
ues converge to unity, as with the smaller system, one in 
particular can be seen to remain much larger than this 
during the entire dynamics. This is a clear indicator of a 
significant SF region within the state as it is dynamically 
driven into the MI regime. For the larger inhomogeneous 
system we see that the slow ramping profile is not adi- 
abatic enough to bring the system into the MI ground 
state shown in Fig. 0{c) . 

Indeed to examine the nature of the state generated by 
the dynamics additional plots of the one-particle density 
matrix are also shown in Fig.0for the two times indicated 
during the dynamics. These show that the state of the 
system remains close in form to that of the ground state 
for the intermediate regime shown in Fig.^Jb), where a 
large SF region exists at the center of the trap. However, 
unlike that ground state we see that sizable correlations 
between the separated SF regions, which were alluded to 
earlier in Sec. IIIII B. do exist and remain present even at 
the bottom of the profile for Fig.^ii) where U/2J = 20. 

Another important difference between the spectra of 
the small and larger system is the more prominent exci- 
tations which have been induced during the transition. 
These are visible as the oscillatory behavior of the eigen- 
values seen in the latter section of the profile in Fig. 
Their presence is consistent with the fact that larger sys- 
tems have more numerous and closely spaced low-lying 
excitations, however, despite this the oscillations have 
only a small amplitude and so do not destroy the SF 
obtained at the end of the transition. 

In order to examine the speed at which coherence is 
reestablished in the system during the latter half (t > t c ) 
of the ramping profile the correlation length of the system 
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FIG. 7: The spectrum of the one-particle density matrix e 7 
for the slow dynamics of the larger inhomogeneous system, 
showing only the largest ten eigenvalues. The dashed lines 
denote the two times (i) t — t c /2 and (ii) t — t c for which 
the absolute value [p m ,n| of the one-particle density matrix is 
plotted. 



must be computed over time. The correlation length is 
typically defined as the distance at which the off-diagonal 
elements of the one-particle density matrix become neg- 
ligible 44] . For symmetrical systems, like those consid- 
ered here, it is natural to measure this from the central 
site m = 25. However, the inhomogeneity of the system, 
which results in the kind correlations between spatially 
separated SF regions just discussed for Fig. [7{i), makes 
the determination of the correlation length ambiguous. 
Instead we choose to examine a cutoff length £ c where 
the spatial correlations with the central site have a spe- 
cific value |p25,25+£ c | = 1/e m 0.37. This value is large 
enough that it corresponds to tracking a point on the cen- 
tral SF region and so can provide a relative measure of its 
size. The change in £ c over time is plotted in Fig. |HJa). 
The same plot also shows the fitted curve whose function 
is that of a smooth 'well', which is the reflection of the 
smooth 'box' used earlier in Eq. H12|) about the t axis. 
The time parameters of this fit are very close to those of 
the resulting 'well' for J(t) generated from the Vo(t) pro- 
file. The variation in £ c over time demonstrates that it is 
capturing the essential changes in the central SF region, 
including the oscillations in its size at later times caused 
by excitations. 

To investigate how much the presence of residual SF 
correlations effects the time dependence of £ c the latter 
half (t > t c ) of the same slow ramping profile was ap- 
plied to the MI ground state of Fig.^fc) found earlier for 
U/2J = 20. In contrast to the state which is dynamically 
driven to the MI regime this ground state has virtually 
no off-diagonal correlations for any site. The change in £ c 
of time for this case is shown in Fig.^b), and as before 
the fitted smooth 'well' function is also plotted. These 
show that despite £ c starting from a much smaller value 
for the MI ground state it still acquires the mean value 
of the final SF state on the same time scale as that of the 
dynamically driven state. 

An important and experimentally motivated measure 
of the coherence of a state can be obtained from its mo- 
mentum distribution function pk- In experiments the in- 
terference pattern resulting from the state of the system 
is examined by allowing all the atoms within the lat- 
tice to expand freely for a short period of time and then 
measure the absorption I{x) at points a; on a distant ob- 
servation line. In the simplest model of this process one 
can neglect both the interactions between atoms during 
the expansion and the spatial dependence of the inter- 
ference caused by the freely evolved Wannier function 
envelopes w(x — x m ) |45j| . This gives the generic fea- 
tures of the interference pattern at an observation point 
x in terms of the path phases acquired by each site in 
a ID lattice of phase coherent matter wave sources. In 
the far-field approximation the intensity I(x) along the 
observation line is proportional to the momentum distri- 
bution p k cx J2 m , n exp[«/c(m - n)]p mi „ HHi^. 

The form of the momentum distribution is sufficiently 
well behaved that its width f w can be determined most 
easily by taking its standard deviation. The variation of 
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FIG. 8: Slow dynamics of the larger inhomogeneous system 
with the variation in the correlation cutoff f c , measured in 
lattice sites, for (a) the complete slow ramping profile starting 
from the SF ground state, and for (b) the MI ground state 
beginning at the bottom of the ramping. The momentum 
distribution width f w , also measured in lattice sites, for the 
same situations is shown in (c) and (d) respectively. The 
dashed curves present in all plots are the fitted smooth 'box' 
or 'well' functions to the data. 



f w in time for the full dynamics is shown in Fig. IHfc), 
along with the fitted smooth 'box' function. The time 
parameters of this fit are again very similar to those of 
J(t). As expected it is seen that f w increases in line with 
the decrease in off-diagonal correlations. The time pro- 
file of f w for the half-ramping of the MI ground state 
is shown in Fig. QJd) and again confirms that the mo- 
mentum distribution width of the SF is reestablished on 
approximately the same time scale as that of the dynam- 
ically driven state. 

Finally we examine the speed at which the correlation 
cutoff length £ c increases with time over the latter half of 
the slow ramping. At any given time t the characteristic 
time scale at which sin gle- atom hopping occurs is given 
by Ttunnci(i) = n/2J(t) |4f|. The simplest description of 
the growth of the central SF region is based on atoms at 
the edge of the system hopping towards the center. In 
this way correlations can be established over the whole 
lattice of M sites Q]. An estimate for the overall time 
scale for this mechanism to occur is given by < rcs toro = 
Afrtunnel/2 0, which for the system and depths used 
here has a value Restore ~ 23 ms. In Fig. |5] the speed 
d^ c /dt obtained from the function fitting is plotted for 



(i) the full dynamics of the ramping and (ii) the half- 
ramping from the MI ground state. In line with the plots 
of £ c in Fig. 0a) and Fig. 0b) we see that there is a time 
delay before there is a significant rate of change in £ c 
for the MI ground state ramping. Since the restoration 
of correlations occurs over the same total time scale in 
both cases the peak in the correlation speed is higher 
for the MI ground state. In addition to these curves the 
characteristic tunnelling speed i>t U nnci(£) = l/itunnei(i) 
over the ramping is also shown, and most importantly 
we note that neither of the two correlation speeds (i) nor 

(ii) exceeds this curve. This confirms that the ramping 
applied is sufficiently slow that the propagation of the SF 
is dominated by single atom hopping. 
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FIG. 9: A comparison between the speeds at which the cor- 
relation cutoff length £ c changes in time for (i) the full slow 
ramping profile starting from the SF ground state and (ii) the 
last half of the ramping profile starting from the MI ground 
state at the bottom of the ramping. The time profile for the 
single-atom hopping speed Vtmnelif) is a l so shown. All speeds 
are expressed in lattice sites per ms. 



B. Fast dynamics 

1. Profile for the fast dynamics 

The time scale over which the slow ramping occurs is 
of the order of 60 ms and so greater than i r ostorc- Here we 
consider ramping occurring much more rapidly. Specifi- 
cally we replace the latter part (t > t c ) of the slow ramp- 
ing profile with a linear ramping of the optical depth 
Vo(t) from Vmi to Vsf as 



V (t) = Vmi 



{Vmi — Vsf) 



(t-tc), 



(13) 



where t runs from t c to t c + t lamp , and the total time of 
the ramping is t ramp . This gives a total ramping profile 
similar to that studied experimentally by Greiner et al 0] . 
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2. Fast dynamics of the larger system: M — 49 

For the fast dynamics we restrict our attention to the 
state that is dynamically driven to the MI regime by the 
slow ramping profile at t = t c . A number of simulation 
runs were performed for total ramping times t ramp be- 
tween 0.1 ms and 10 ms, along with t ramp = ms which 
is equivalent to the initial state. The value of £ c ob- 
tained at the end of each of the ramping times is plotted 
in Fig. HOf a) . We see that there is a steady monotonic 
increase in the £ c for ramping times i ramp , except where 
it is broken by peaks and troughs which are the expected 
manifestations of the trapping used. In particular the 
trough centered around t lamp ~ 7 ms corresponds to the 
period of an oscillation with frequency 2uj, where u> is the 
trapping frequency introduced in Sec. lIIII B. On a similar 
basis the spikes which appear around t ramp ~ 1 rns can 
be seen to be a result of the excitation spectrum. 
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FIG. 10: Results for rapid dynamics, (a) the final correlation 
cutoff length £ c obtained for different linear ramping times 
tramp, focusing on f ra mp = 8 ms we have (b) the spectrum 
of the one-particle density matrix e 7 showing only the largest 
ten eigenvalues, and (c) the variation in £ c over the simulation 
run, along with a fitted smooth 'box' function shown as the 
dashed curve. 



We take a special interest in the ramping time t ramp = 
8 ms where £ c obtains its maximum value approximately 
equal to that of the SF ground state. The variation in 



time of the spectrum e 7 of the one-particle density matrix 
and £ c during this particular ramping simulation is given 
in Fig. llOf b) and Fig. HOf c). A well-behaved monotonic 
increase in £ c is observed which can be accurately fitted 
over the interval [0,8] ms by a smooth 'box' function, 
as used earlier. This again provides the basis for com- 
puting the speed d^ c /dt at which the correlation cutoff 
length £ c is increasing over the ramping and is shown in 
Fig. lllf b) along with that of the characteristic tunnelling 
speed ftunnei(i) for the rapid ramping profile. Unlike the 
similar comparison for the slow dynamics we see here that 
after approximately 3 ms £ c is increasing in time much 
more rapidly than the single-site tunneling speed i>t U nnei 
alone can account for. Indeed by the end of the ramping 
d^ c /dt is almost three times that of the maximum tun- 
neling speed. This is a clear indication that single atom 
hopping is not adequate to describe the growth of the 
central SF region for such rapid dynamics. Instead the 
specific form and contributions of higher order correla- 
tion functions must play a crucial role. 
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FIG. 11: Rapid ramping of an optical lattice from the MI 
(U/2J = 20) to the SF (U/2J = 2) regime for N = 40 atoms 
in M = 49 lattice sites superimposed by a magnetic trapping 
potential. The width of the central interference fringe f w as a 
function of the ramping time tramp is shown in (a). The solid 
curve is a fit using a double-exponential decay (n = 0.22 ms, 
T2 — 2.14 ms) (cf. In (b) the rate of change of the cor- 
relation cutoff length £ c is shown for the ramping performed 
with tramp = 8 ms, along with the profile for the characteristic 
tunneling speed vtunael(t) for the ramping. Both are plotted 
in units of lattice sites per ms. 



To draw direct comparisons with the results of Grcincr 
et al. on the restoration of coherence we also plot the 
momentum distribution width f w obtained from the one- 
particle density matrix at the end of each ramping in 
Fig. lllf ah The data points for this quantity show a pro- 
nounced trend, without any of the large-scale variations 
seen in £ c caused by the trapping |4l|. The decrease of 
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f w with increasing i ramp fits well to a double-exponential 
decay curve of the form 

f w (Wp) = A 1 e -*Wn + A 2 e -WM + c> (14) 

where t%, t 2 are the characteristic decay times, A\, A 2 are 
coefficients, and C is a constant. Most notably the exact 
same functional form was found to fit the experimental 
data in 1]. Since their experiment was conducted for 
a 3D lattice, along with a larger mean occupancy and 
a deeper ramping profile, exact agreement for the time 
parameters of this fit is not expected. However, we do 
note that the ratio of the decay times used for their fit 
and ours are both t 2 /ti ss 10. Similarly we can make the 
same observation as made in Q that the momentum dis- 
tribution width f w has returned to its steady state value 
within a time scale approximately of order Ttunnel- This 
is much shorter than the expected time trcstoro required 
for coherence to spread over the whole lattice of M sites 
via single-atom hopping. This confirms that the restora- 
tion of coherence as seen in the experiment is accurately 
described by the BHM. 



3. Validity of the simulation for fast dynamics 



is reestablished within the system for both slow and rapid 
dynamics which cross the SF-M1 transition. Our results 
indicate that for slow ramping of the lattice depth the SF 
growth is consistent with single atom hopping as might 
naively be expected. However, for very rapid ramping 
of the lattice depth we find that the SF growth is much 
greater than can be explained by this mechanism alone 
and so points to the importance of higher-order correla- 
tion functions. We made direct comparisons between our 
simulation results for the momentum distribution width 
f w during rapid ramping and the experimental results 
obtained by Greiner et al. [l| and found that the reduc- 
tion in f w with the ramping time follows precisely the 
same functional form as their data, despite a number 
of significant differences in the systems analyzed. Per- 
haps most fundamentally we have shown that the re- 
sults obtained in for the rapid restoration of coher- 
ence are consistent and explicable within the BHM alone 
and are present even in ID systems. Finally, we note 
that a detailed knowledge of the correlations of atoms 
in different sites and the particle number fluctuations 
as provided by our numerical calculations are impor- 
tant for utilizing the MI state in a number of applica- 

tions nasiiiiammEH. 



The simulations performed here assume that the dy- 
namics of the atoms is described by the lowest Bloch 
band of the optical lattice. This assumption holds if the 
typical frequency / ss l/t ramp of the ramp in U and J 
obeys where v = ^JAErVq/2'k is the harmonic 

approximation of the excitation frequency to the first ex- 
cited Bloch band |15j . The shortest ramping time we 
considered is t tamp = 0.1 ms, while \ jv = 0.05 ms for 
the lattice on average over the ramping. Because the 
condition / <C v is not fulfilled, we numerically calcu- 
lated the probability of exciting a single particle, when 
initially prepared in the lowest Bloch band and located at 
the central site of the lattice, during the ramping above 
as a function of t ramp . We find that the time evolution 
is well approximated by the adiabatic time evolution for 
^ramp > 0.05 ms and that it changes to being sudden for 
iramp < 0.005 ms. Therefore we expect only a small in- 
fluence to the form of the curve in Fig. Illf a) between 
the points at t ram p = ms and t ramp = 0.1 ms due to 
higher band excitations which is not resolved in the ex- 
periments [l|. 



V. CONCLUSION 

In these studies we have established the accuracy and 
applicability of the TEBD algorithm to the BHM, for 
both the computation of ground states and its dynam- 
ics. We have then applied this method to systems of a 
size equivalent to those studied in experiments and in 
the presence of a trapping potential. In particular we 
have examined the nature and speed in which coherence 
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APPENDIX A: ID OPTICAL LATTICES 

The starting point for our physical model is the Hamil- 
tonian for weakly interacting bosonic atoms in an exter- 
nal trapping potential 



II I drV> f (r 
1 iira 



1 



2771,4 



V 2 + V (r) + V t (t) ^(r) 



2 rriA 



drip' (r)V^ (r)ip(r)ip(r) , 



(Al) 



with ip (r) the bosonic field operator for atoms in a given 
internal atomic state, Vo{r) is the optical lattice poten- 
tial, and Vr(j) describes a slowly varying additional ex- 
ternal trapping potential such as that created by mag- 
netic fields. The interaction between the atoms is mod- 
eled by a contact potential with s-wave scattering length 
a s and ttia is the mass of the atoms. 

We assume the optical lattice potential to have the 
form Vo(r) = X)j=i ^S'o sin 2 (/crj) with wave number k = 
2ir/\ and A the wavelength of the laser light yielding a 
lattice period a = A/2. The spatial coordinate is de- 
noted by r = (ri,r 2 ,r 3 ). This lattice potential can be 
realized by interfering three pairs of counterpropagating 
laser beams from three orthogonal directions. The height 
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of the potential Vjo is proportional to the intensity of the 
lasers in the j th pair of laser beams. We assume the 
intensity to be very large in the ri and r% directions so 
that the atoms do not tunnel in either of these directions. 
Hence their motion is restricted to the r% = x direction 
and this optical lattice setup thus allows the creation 
of effective ID systems. The resulting atomic pipelines 
0- EH H3 are we ll isolated from each other and we can 
thus restrict our considerations to just one of them. 

The center position of the lattice site to of this ID 
system is given by x rn — ma, and so a particle occu- 
pying the lowest Bloch band which is localized at this 
site is described by the wave function <j) m (r) = w\ (x — 
x m )w2{T2)wz{rz), where Wj are the Wannier functions of 
the lowest Bloch band in the j th direction. By ne- 
glecting all excitations to higher bands and expanding the 
bosonic field operators into the mode functions <j) m (r) the 
Hamiltonian H reduces to the ID Bose-Hubbard model 
[l5| given in Eq. (JIJ Sec.[H]A . The parameter U of the 
BHM is given by U — Aira s J dr|0 m (r)| 4 /itia and cor- 
responds to the strength of the on-site repulsion of two 
atoms occupying the same lattice site to. The hopping 
matrix element J between adjacent sites m and m + 1 is 
given by 

+ V sm 2 (kx)^wi(x - x m+ i). (A2) 

The numerical values for U and J for different depths of 
the optical lattice Vo = Via can be found in |T^ . 

APPENDIX B: MF APPROXIMATION AND 
GUTZWILLER ANSATZ 

The presence of a trapping potential causes the sys- 
tem to become inhomogeneous allowing the coexistence 



of spatially separated SF and MI regions E IHIH |H 
Eol [4ll lia |43| . Now each lattice site has a local chem- 
ical potential fi m , and together with the overall U/2J 
parameter for the system this gives a phase diagram co- 
ordinate p = (2 J/U, Hm/U). The MF determination of 
which regime a given lattice site lies in is then based on 
where preci sely p resides in the homogeneous BHM phase 
diagram |5l l4ll |48| . 

The MF calculations performed utilized the standard 
decoupling of the hopping term 0, 



btfim+i = b\ n (b m+1 ) + (b1 n ) b m+1 - (hi) (b m+1 ) , (Bl) 



which decompose the BHM Hamiltonian in Eq. into 
a set of M single site Hamiltonians H m . Specifically the 
single site H m has the form 



H m = -2J((j) m bl n +h.c.)-2\(j) m \ 2 -^ m bl l b m +—bl n bl n b m b m , 

(B2) 

for each lattice site to, dependent on the sites superfluid 
order parameter cj> m = (b m ), which is assumed to vary 
slowly over the system, and its local chemical potential 
fj, m . The MF ground state is then a product state over 
sites |^mf) = Ilm=i l^™) determined by minimizing the 
complete set of M-site Hamiltonians H m with respects 
to the order parameters (f> m and then diagonalizing to 
extract eigenstate \ip™} with the smallest eigenvalue for 
each site. Given that each [ip™) = ^2 n c„ m \n m ) this 
procedure is equivalent to approximation methods based 
on the Gutzwiller ansatz 0, and has been applied 
successfully in modeling the qualitative behavior of BHM 
phase diagram |48| . 
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